Modularbeit Wiener Filter¶
Herleitung des Wiener Filters:¶

Suche: $$\begin{equation*} \underset{W}{argmin} \; \mathbb{E}\Big[\big\|\hat f - f\big\|^2\Big] \end{equation*}$$
Annahmen:¶
$$\begin{align*} g = h \ast f + n \; &\Longrightarrow \; G = H \cdot F + N \tag{1}\\ \hat f = g \ast w \; &\Longrightarrow \; \hat F = G \cdot W \tag{2}\\ \end{align*}$$
$$\begin{align*} S_{FF} &= \mathbb{E}\Big[\big\|F\big\|^2\Big] \tag{3}\\ S_{NN} &= \mathbb{E}\Big[\big\|N\big\|^2\Big] \tag{4}\\ \end{align*}$$
Das Bild und das Rauschen sind unkorreliert: $$\begin{align*} S_{FN} = \mathbb{E}\Big[F^*N\Big] = \mathbb{E}\Big[FN^*\Big] = 0 \tag{5}\\ \end{align*}$$
Herleitung:¶
Minimierung des quadratischen Fehlers vom Orginalbild und des rekonstruierten Bildes: $$\begin{equation*} \mathbb{E}\Big[\big\|\hat F - F\big\|^2\Big] \overset{(2)}{=} E\Big[\big\|GW - F\big\|^2\Big] \\ \end{equation*}$$
$$\begin{equation*} = \mathbb{E}\Big[\big\|GW\big\|^2 + \big\|F\big\|^2 - (GW)^*F - GWF^*\Big] \\ \end{equation*}$$
$$\begin{equation*} = \big\|W\big\|^2 \,\mathbb{E}\Big[\big\|G\big\|^2\Big] + \mathbb{E}\Big[\big\|F\big\|^2\Big] - W^* \,\mathbb{E}\Big[G^*F\Big] - W \,\mathbb{E}\Big[GF^*\Big] \\ \end{equation*}$$
$$\begin{equation*} \overset{(1,3)}{=} \big\|W\big\|^2 \,\mathbb{E}\Big[\big\|HF + N\big\|^2\Big] + S_{FF} - W^* \,\mathbb{E}\Big[(HF + N)^*F\Big] - W \,\mathbb{E}\Big[(HF + N)F^*\Big] \\ \end{equation*}$$
$$\begin{equation*} = \big\|W\big\|^2 \,\mathbb{E}\Big[\big\|HF\big\|^2 + \big\|N\big\|^2 + (HF)^*N + HFN^*\Big] + S_{FF} - W^* \,\mathbb{E}\Big[H^*\big\|F\big\|^2 + N^*F\Big] - W \,\mathbb{E}\Big[H\big\|F\big\|^2 + NF^*\Big] \end{equation*}$$
$$\begin{align*} = &\big\|W\big\|^2 \Big(\big\|H\big\|^2 \,\mathbb{E}\Big[\big\|F\big\|^2\Big] + \mathbb{E}\Big[\big\|N\big\|^2\Big] + H^* \,\mathbb{E}\Big[F^*N\Big] + H \,\mathbb{E}\Big[FN^*\Big]\Big) + S_{FF} \\ &- W^* \Big(H^* \,\mathbb{E}\Big[\|F\|^2\Big] + \mathbb{E}\Big[N^*F\Big]\Big) - W \Big(H \,\mathbb{E}\Big[\big\|F\big\|^2\Big] + \mathbb{E}\Big[NF^*\Big]\Big) \\ \end{align*}$$
$$\begin{equation*} \overset{(3,4,5)}{=} \big\|W\big\|^2 \Big(\big\|H\big\|^2 S_{FF} + S_{NN} + H^*0 + H0\Big) + S_{FF} - W^* \Big(H^* S_{FF} + 0\Big) - W \Big(H S_{FF} + 0\Big) \\ \end{equation*}$$
$$\begin{equation*} = \big\|W\big\|^2 \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) + S_{FF} - W^*H^*S_{FF} - WHS_{FF} \end{equation*}$$
Ableiten und Nullstellen: $$\begin{align*} \frac{\partial \, \mathbb{E}\Big[\big\|\hat F - F\big\|^2\Big]}{\partial \, W} &= \frac{\partial \, \big\|W\big\|^2 \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) + S_{FF} - W^*H^*S_{FF} - WHS_{FF}}{\partial \, W} \\ &= W^* \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) - HS_{FF} \overset{!}{=} 0 \\ \end{align*}$$
$$\begin{align*} &\Longleftrightarrow W^* \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) = HS_{FF} \\ &\Longleftrightarrow W \Big(\big\|H\big\|^2 S_{FF} + S_{NN}\Big) = H^*S_{FF} \\ &\Longleftrightarrow W = \frac{H^*S_{FF}}{\|H\|^2 S_{FF} + S_{NN}} = \frac{H^*}{\|H\|^2 + \frac{S_{NN}}{S_{FF}}} \end{align*}$$
Implementation:¶
from PIL import Image
import numpy as np
from matplotlib import pyplot as plt
from tqdm import tqdm
import math
def read_image(image_path):
image = np.array(Image.open(image_path).convert('L'))
normalized_image = image / image.max()
return normalized_image
def show_image(image):
plt.imshow(image, cmap='gray')
plt.axis('off')
plt.show()
def fft(image):
return np.fft.fft2(image)
def ifft(image_fq):
return np.abs(np.fft.ifft2(image_fq))
Mean Squared Error:¶
$$\begin{equation*} MSE = \mathbb{E}\Big[\big\| \hat f - f \big\|^2\Big] \end{equation*}$$
Peak Signal To Noise Ratio:¶
$$\begin{equation*} PSNR = 10\log_{10}\Bigg(\frac{MAX_f^2}{MSE}\Bigg) \end{equation*}$$
Accuracy Amoothness Error:¶
$$\begin{equation*} ASE = \mathbb{E}\Big[\big\| \hat g - h \ast \hat f \big\|^2\Big] + \lambda \, \mathbb{E}\Big[\big\| \hat f^{\prime} \big\|^2\Big] \end{equation*}$$
def mse(f, f_hat):
return (np.square(f - f_hat)).mean()
def psnr(f, f_hat):
mse_value = mse(f, f_hat)
if mse_value == 0:
return 0
return 10 * np.log10(1 / mse_value)
def accuracy_smothness_error(g, f_hat, H, _lambda):
data_accuracy = np.linalg.norm(g - ifft(fft(f_hat) * H))
smothness = np.sum(np.square(np.gradient(f_hat)))
return data_accuracy + _lambda * smothness
Wiener Filter:¶
$$\begin{align*} \hat f = \mathcal{F}^{-1} \Bigg\{\Bigg[\frac{H^*}{\|H\|^2 + K}\Bigg] G \Bigg\} \end{align*}$$
def wiener_filter(G, H, K):
W = np.conjugate(H) / (np.abs(np.square(H)) + K)
F_hat = W * G
f_hat = ifft(F_hat)
return f_hat
Motion Blur Filter:¶
$$\begin{equation*} H(u,v) = T sinc\Big(\pi\big(u \Delta x + v \Delta y\big)\Big)\mathrm{e}^{\mathrm{i} \pi \big(u \Delta x + v \Delta y\big)} \end{equation*}$$
def get_motion_blur_filter_fq(shape, T, dx, dy):
U_DIM ,V_DIM = shape
H_u = np.arange(0, U_DIM, 1)
H_v = np.arange(0, V_DIM, 1)
H_U, H_V = np.meshgrid(H_v, H_u)
H = np.dstack((H_U, H_V))
dxy = np.array([dx, dy])
def s(uv):
return np.pi * uv.dot(dxy)
H = np.apply_along_axis(s, axis=2, arr=H)
H = T * np.sinc(H) * np.exp(-1j * H)
return H
Gausian Noise Filter:¶
$$\begin{align*} n(x,y) &\sim \mathcal{N}\big(\mu,\sigma^2\big)\\ N(u,v) &= \mathcal{F}\Big\{n(x,y)\Big\} \\ \end{align*}$$
def get_gausian_noise_filter_fq(shape, variance, mean=0):
n = np.random.normal(mean, variance, shape)
N = fft(n)
return N
Plotting:¶
def plot_images_restauration(original_image, disturbed_image, restored_images, restored_titels):
fig, ax = plt.subplots(2, math.ceil(len(restored_images) / 2) + 1)
fig.set_size_inches((14, 12))
ax[0,0].set_title(f'Original Image\n MSE: {mse(original_image, original_image):.2f}')
ax[0,0].imshow(original_image, cmap='gray')
ax[0,0].axis('off')
ax[1,0].set_title(f'Disturbed Image\n MSE: {mse(original_image, disturbed_image):.2f}')
ax[1,0].imshow(disturbed_image, cmap='gray')
ax[1,0].axis('off')
for i in range(len(restored_images)):
x = int(i / 2) + 1
y = i % 2
ax[y, x].set_title(f'{restored_titels[i]}\n MSE: {mse(original_image, restored_images[i]):.2f}')
ax[y, x].imshow(restored_images[i], cmap='gray')
ax[y, x].axis('off')
plt.tight_layout()
plt.show()
def optimize_k(disturbed_image_fq, motion_blur_filter_fq, _lambda=0, _dropout=True, _max_iter=50):
disturbed_image = ifft(disturbed_image_fq)
best_error = np.inf
best_k = 0.5
step = 0.5
for _ in tqdm(range(_max_iter), desc='Optimizing K'):
last_iter_error = best_error
K = np.linspace(best_k + step, best_k - step, 5, endpoint=False)[1:]
for k in K:
restored_image = wiener_filter(disturbed_image_fq, motion_blur_filter_fq, k)
error = accuracy_smothness_error(disturbed_image, restored_image, motion_blur_filter_fq, _lambda)
if error < best_error:
best_error = error
best_k = k
step /= 2
if _dropout and last_iter_error - best_error < 0.01:
break
return best_k
Wiener filtering By Averaging Over Noise:¶
[Yoo, Jae‐Chern, and Chang Wook Ahn. "Image restoration by blind‐Wiener filter." IET Image Processing 8.12 (2014): 815-823.]

def wiener_filter_averaging(G, H, variance, n=10):
f_hat = 0
for _ in range(n):
N_m = get_gausian_noise_filter_fq(G.shape, variance)
D_hat_m = G - N_m
S_D_hat_m= np.mean(np.abs(np.square(D_hat_m)))
S_N_m = np.mean(np.abs(np.square(N_m)))
f_hat_m = wiener_filter(D_hat_m, H, S_N_m/S_D_hat_m)
f_hat += f_hat_m
return f_hat / n
def test_wiener_filter_averaging(image_path, _mb_T=1, _mb_dx=0.005, _mb_dy=0.001, _g_var=0.05):
f = read_image(image_path)
F = fft(f)
H = get_motion_blur_filter_fq(F.shape, _mb_T, _mb_dx, _mb_dy)
N = get_gausian_noise_filter_fq(F.shape, _g_var)
G = H * F + N
g = ifft(G)
images = []
titles = []
variances = [0.01, 0.1, 0.2, 0.4, 0.6, 0.8]
for var in tqdm(variances):
images.append(wiener_filter_averaging(G, H, var))
titles.append(f'Averaging (var): {var}')
plot_images_restauration(f, g, images, titles)
test_wiener_filter_averaging('images/schach.jpg', _mb_T=10, _mb_dx=0.01, _mb_dy=0.1, _g_var=0.1)
100%|██████████| 6/6 [00:13<00:00, 2.18s/it]
Other Posibilities Of Blind Wiener Filtering:¶
- Parameter Estimation of Power Spectral Density of Noise and Image
- location depentently K-parametric wiener filtering: K as adaptively varying-coefficient
- ...
Testing:¶
def test_restaurations(image_path, _mb_T=1, _mb_dx=0.005, _mb_dy=0.001, _n_var=0.05):
f = read_image(image_path)
F = fft(f)
H = get_motion_blur_filter_fq(F.shape, _mb_T, _mb_dx, _mb_dy)
N = get_gausian_noise_filter_fq(F.shape, _n_var)
G = H * F + N
g = ifft(G)
# wiener filter with knowledge of powerspectrum of image and noise
S_F = np.mean(np.abs(np.square(F)))
S_N = np.mean(np.abs(np.square(N)))
f_hat = wiener_filter(G, H, S_N/S_F)
# wiener filter with optimization of K
optimal_k = optimize_k(G, H, _lambda=0)
f_hat_2 = wiener_filter(G, H, optimal_k)
# wiener filter with averaging over noise
f_hat_3 = wiener_filter_averaging(G, H, 0.01)
f_hat_4 = wiener_filter_averaging(G, H, 0.2)
plot_images_restauration(f, g, [f_hat, f_hat_2, f_hat_3, f_hat_4], ['Wiener Filter', 'Optimization', 'Averaging (var: 0.01)', 'Averaging (var: 0.2)'])
Only Motion Blur:¶
test_restaurations('images/lion.jpg', _mb_T=1, _mb_dx=0.01, _mb_dy=0.01, _n_var=0)
Optimizing K: 56%|█████▌ | 28/50 [00:03<00:02, 8.81it/s]
Only Noise:¶
test_restaurations('images/schach.jpg', _mb_T=1, _mb_dx=0, _mb_dy=0, _n_var=0.3)
Optimizing K: 26%|██▌ | 13/50 [00:17<00:48, 1.32s/it]
Motion Blur and Noise:¶
test_restaurations('images/schach.jpg', _mb_T=1, _mb_dx=-0.001, _mb_dy=-0.02, _n_var=0.1)
Optimizing K: 6%|▌ | 3/50 [00:05<01:18, 1.67s/it]
test_restaurations('images/eule.png', _mb_T=1, _mb_dx=0.001, _mb_dy=0.001, _n_var=0.2)
Optimizing K: 4%|▍ | 2/50 [00:07<02:56, 3.69s/it]
test_restaurations('images/eule.png', _mb_T=10, _mb_dx=0.001, _mb_dy=-0.01, _n_var=0.01)
Optimizing K: 16%|█▌ | 8/50 [00:17<01:33, 2.23s/it]
test_restaurations('images/schach.jpg', _mb_T=10, _mb_dx=0.1, _mb_dy=0.2, _n_var=0.1)
Optimizing K: 6%|▌ | 3/50 [00:04<01:14, 1.59s/it]